
#!/usr/bin/perl
use strict;
use Statistics::Distributions; # For ChiSquare P-value calculation
use Text::NSP::Measures::2D::Fisher2::twotailed; # For Fisher Exact Test P-value

my $tpeddir = "/data1/bsi/refdata/genetics/1000genomes/processed_data/pilot_data/2010_03/pilot_1/CEU/no_triallelic/";
my $tgdir = "/data1/bsi/refdata/genetics/1000Genomes/ftp/pilot_data/release/2010_03/pilot1/";
my $tgfile = "CEU.SRP000031.2010_03.genotypes.vcf.gz";

#my @chrs = ("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22");

# Define trustable LD region. (e.g. how far away to look).
my $range = 1000000; # 1 Megabytes

# Define target region
my $region_start = 61741000;
my $region_end = 61858000;
my $region_chr = 20;
#my $region_start = 130400000;
#my $region_end =   130800000;
#my $region_chr = 8;

my @region_snps =();
my @region_snps_pos =();
my @target_allele =();
my @target_allele_type =();
my @target_allele_id=();


# INSERT SOME CODE TO FILL UP the region_snps array from the region information
# Can either can the tped files at $tpeddir or read the SNPs from an array annotation file.
my $tpeddir = "/data1/bsi/refdata/genetics/1000Genomes/processed_data/pilot_data/2010_03/pilot_1/CEU/no_triallelic/";
#22 rs62224609 0 14431249 T T T T T T T T T T T T T T T T T C T T T T T T T T T T T T T T T T T T T T T C T T T T T T T T T T T T T T T C T T T T T T T C T T T T T T T T T T T T T T T T T T T T T T T C C T T T T C T T T T T T T T T T T C T C T T T T T T T T T T T T 

push(@region_snps,"rs6010620");
push(@region_snps,"rs2297440");
push(@region_snps,"rs4809324");

#push(@region_snps,"rs10464870");
#push(@region_snps,"rs891835");
#push(@region_snps,"rs6470745");
#push(@region_snps,"rs16904140");
#push(@region_snps,"rs4295627");

push(@region_snps_pos,61780283);
push(@region_snps_pos,61782743);
push(@region_snps_pos,61788664);
#push(@region_snps_pos,130547005);
#push(@region_snps_pos,130560934);
#push(@region_snps_pos,130711103);
#push(@region_snps_pos,130734825);
#push(@region_snps_pos,130754639);

my $tpos = ((sort {$a<=> $b} @region_snps_pos)[0])  . "-" .   ((sort {$a<=>$b} @region_snps_pos)[$#region_snps_pos]);
if($#region_snps_pos==0) {
    $tpos = $region_snps_pos[0];
}

# Target Allele on Plus strand
#push(@target_allele,'');
#push(@target_allele,'');
push(@target_allele,'G');
push(@target_allele,'C');
push(@target_allele,'T');

# The alleles are coded as 0 (reference) and 1 (alternate)
for(my $i=0;$i<=$#target_allele;$i++) {
    push(@target_allele_id,0); # reference
}
# either L/H (Low or High MAF) or R/A (reference/Alternate)
#push(@target_allele_type,'L'); 
#push(@target_allele_type,'L');
push(@target_allele_type,'A');
push(@target_allele_type,'R');
push(@target_allele_type,'R');

# The alleles are coded as 0 (reference) and 1 (alternate)
for(my $i=0;$i<=$#target_allele;$i++) {
    if($target_allele_type[$i] eq 'R') {
	push(@target_allele_id,0); # reference
    } elsif($target_allele_type[$i] eq 'A') {
	push(@target_allele_id,1); # Alternate
    } else {
	push(@target_allele_id,-1); # unknown 
    }
}


my $nregion_snps = 1 + $#region_snps;



my %REFS; # Store the vector of Haplophase-specific 1's if the haplotype == the reference allele
my %ALTS;  # Store the vector of Haplophase-specific 1's if the haplotype == the alternate allele
my %INFO; # $INFO{$snpid} = \($ac,$an,$na,$data[$iref],$data[$ialt],$maf); 
my %GENO; # For each SNP, stores a ref to a list of genotypes(The count of the major allele)
my %REFAL;
my %ALTAL;
my %HAPL;  # Haplotype flag Y/N/? to indicate if the target variant contains the target allele on each haplotype.


# Open VCF file and scan down file.
# Load 


open FIN, "gunzip -c $tgdir$tgfile | ";

my @head;
my $n0=0;
my $l;
my $iinfo=-1; # indec for info
my $iref=-1; #index for ref allele - genotype code ==0
my $ialt=-1; # index for alt allele - genotype code ==1 
my $ifirsts=-1; # index of first sample
my $ichr=-1; # Should be 0
my $ipos=-1; # Should be 1
my $isnpid=-1;  # Should be 2 if any
my $nfound=0;

my $nchroms=0;

while(($l=<FIN>)  && ($nfound<$nregion_snps)) {
    $l =~ s/[\r,\n,\l]+$//g;
    if($l =~ /^#.*/) {
       if($l =~ /^#CHROM.*/){
	  $l =~ s/^#//;
	  @head = split(/\t/,$l);
	  for(my $i=0;$i<=$#head && $ifirsts==-1; $i++) {
	      if($head[$i] =~ /^INFO$/) {
		  $iinfo=$i;
	      } elsif($head[$i] =~ /^REF$/) {
		  $iref=$i;
	      } elsif($head[$i] =~ /^ALT$/) {
		  $ialt=$i;
	      } elsif($head[$i] =~ /^ID$/) {
		  $isnpid=$i;
	      } elsif($head[$i] =~ /^CHROM$/) {
		  $ichr=$i;
	      } elsif($head[$i] =~ /^POS$/) {
		  $ipos=$i;
	      } elsif($ifirsts==-1 && $head[$i] =~ /^NA.*/) {
		  $ifirsts=$i;
	      }
	  }
      }
   } else {
	    my @data = split(/\t/,$l);
	    my $snpid = $data[$isnpid];
	    my $foundit=-1;
	    for(my $j=0;$foundit==-1 && $j<=$#region_snps;$j++) {
		if($region_snps[$j] eq $snpid) {
		    $foundit=1;
		    $nfound++;
		    my $f=$data[$iinfo];
	      
		    my @info = split(";",$f);

		    my $ac=0;
		    my $an=0;
		    for(my $i=0;$i<=$#info;$i++) {
			my @kv = split("=",$info[$i]);
			if($kv[0] =~ /AC/) {# Non-reference Allele Count
						$ac=$kv[1];
					    }elsif($kv[0] =~ /AN/) { # Total Chromosome Count
						$an=$kv[1];
					    }
		    }
		    
		    my $maf=0;
		    if($an>0) {
			if($ac> ($an-$ac)) {
			    $maf=($an-$ac)/$an; # Reference Allele is minor Allele
			    if($target_allele_type[$j] eq 'L') {
				$target_allele_id[$j]=0; # alternate
			    } elsif($target_allele_type[$j] eq 'H') {
				$target_allele_id[$j]=1; # Reference is common Alleles.
			    }
			} else {
			    $maf=$ac/$an; # Alternate Allele is minor Allele
			    if($target_allele_type[$j] eq 'L') {
				$target_allele_id[$j]=1; # alternate
			    } elsif($target_allele_type[$j] eq 'H') {
				$target_allele_id[$j]=0; # Reference is common Alleles.
			    }
			    
			}
		    }
		    my @refallele=();
		    my @altallele=();
		    my @genos=();
		    my $na=0; # Number of ambigous;
		    my @has_target_allele=(); # 'Y'  if the haplotype contains the target allele. N if not,'?' is not sure.
		    for(my $i=$ifirsts;$i<=$#data;$i++) {
			my @genoarr = split(":",$data[$i]);
			my $genostring = $genoarr[0];
		    # SAVE INFORMATION
			if($genostring =~ /^([0-2])\|([0-2])$/) {
			    my $al1 =$1;
			    my $al2 =$2;
			    push(@genos,$al1+$al2);
			    if($al1 ==0) {
				push(@refallele,1);
				push(@altallele,0);
			    } else {
				push(@refallele,0);
				push(@altallele,1);
			    }
			    if($al1 == $target_allele_id[$j]) {
				$has_target_allele[2*($i-$ifirsts)]='Y';
			    } else {
				$has_target_allele[2*($i-$ifirsts)]='N';
			    }
			    if($al2 ==0) {
				push(@refallele,1);
				push(@altallele,0);
			    } else {
				push(@refallele,0);
				push(@altallele,1);
			    }
			    if($al2 == $target_allele_id[$j]) {
				$has_target_allele[2*($i-$ifirsts)+1]='Y';
			    } else {
				$has_target_allele[2*($i-$ifirsts)+1]='N';
			    }
			} else { # Ambigous phase (e.g. genotype called as "/"
			    push(@refallele,0); # When both are missing, treated as missing data
			    push(@altallele,0);
			    $na++;
			    if($genostring =~ /^([0-2])\/([0-2])$/) {
				push(@genos,$1+$2);
				my $al1 =$1;
				my $al2 =$2;
				if(($al1 != $target_allele_id[$j]) && ($al2 != $target_allele_id[$j])) {
				    $has_target_allele[2*($i-$ifirsts)]='N';
				    $has_target_allele[2*($i-$ifirsts)+1]='N';
				} else {
				    $has_target_allele[2*($i-$ifirsts)]='?';
				    $has_target_allele[2*($i-$ifirsts)+1]='?';
				}
				    
			    } else {
				push(@genos,-1); # Missing data
				$has_target_allele[2*($i-$ifirsts)]='?';
				$has_target_allele[2*($i-$ifirsts)+1]='?';
			    }

			}
		    }
		    $HAPL{$snpid}=\@has_target_allele;
		    $nchroms = 1+$#has_target_allele;
		    $REFS{$snpid}=\@refallele;
		    $ALTS{$snpid}=\@altallele;
		    $GENO{$snpid}=\@genos;
		    my @snpinfo = ($ac,$an,$na,$data[$iref],$data[$ialt],$maf,$data[$ipos],$data[$ichr]); 

		    $INFO{$snpid} = \@snpinfo;
		    $REFAL{$snpid} = $data[$iref];
		    $ALTAL{$snpid} = $data[$ialt];
		}
	    }
	}
}
close FIN;
my @hasHaplotype;
my @foundSNPs = keys %HAPL;
for(my $i=0;$i<$nchroms;$i++) {
    push(@hasHaplotype,-2); # By default, -1 means unknown.
}
for(my $j=0;$j<=$#foundSNPs;$j++) {
    my $hta_ref = $HAPL{$foundSNPs[$j]};
    my @hta = @$hta_ref;
    for(my $i=0;$i<$nchroms;$i++) {
	if($hta[$i] eq 'N') { # A single No is enough to make this NOT the haplotype
	    $hasHaplotype[$i]=0;
	} elsif(($hta[$i] eq 'Y') && $hasHaplotype[$i]==-2) {
	    $hasHaplotype[$i]=1;
	} elsif(($hta[$i] eq '?') && $hasHaplotype[$i]==1) {
	    $hasHaplotype[$i]=-1;
	}
    }
}
for(my $i=0;$i<$nchroms;$i++) {
    if($hasHaplotype[$i]==-2) {
	$hasHaplotype[$i]=-1;
    }
}
my $indiv_with_hapl=0;
my %INDIVHASHAPL;
for(my $i=0;$i<$nchroms;$i++) {
    if($hasHaplotype[$i]==1) {
	$indiv_with_hapl++;
	my $indiv = $head[$ifirsts+($i>>1)];
	if(exists $INDIVHASHAPL{$indiv}) {
	    $INDIVHASHAPL{$indiv}=2;
	} else {
	    $INDIVHASHAPL{$indiv}=$i-2*($i>>1);
	}
    }
}

my $tmaf = $indiv_with_hapl/(1+$#hasHaplotype);
die "Did not find any INDIV with target haplotype" unless $indiv_with_hapl>0;
#{
#print "The following inviduals had the target haplotype ";
#my @samples = keys %INDIVHASHAPL;
#for(my $i=0;$i<=$#samples;$i++) {
#    my $hapcode = $INDIVHASHAPL{$samples[$i]};
#    if($hapcode ==2) { print $samples[$i] . " has two copies\n";}
#    else { print $samples[$i] . " has 1 copy\n";}
#}

#}
#die "";
#print "# Found $nfound of the " . (1+$#region_snps). " target SNPs\n";
#my %REFS; # Store the vector of Haplophase-specific 1's if the haplotype == the reference allele
#my %ALTS;  # Store the vector of Haplophase-specific 1's if the haplotype == the alternate allele
#my %INFO; # $INFO{$snpid} =\($ac,$an,$na,$data[$iref],$data[$ialt],$maf,$data[$ipos],$data[$ichr]); 

print "target_snpid\ttarget_maf\tchr\ttarget_pos\ttref_allele\ttalt_allele\tsnpid\tmaf\tpos\tref_allele\talt_allele\tr2_ld\tr2_pearson\tchisq\tchisq_pvalue\tfisher_two_tailed\tfisher_missing_in_test\tfisher_missing_in_target\tfisher_missing_causes_missing\tn_tref_ref\tn_tref_alt\tn_talt_ref\tn_talt_alt\tn_tref_na\tn_talt_na\tn_na_ref\tn_na_alt\tn_na_na\n";


my $first_pos = $region_start - $range;
my $last_pos = $region_end + $range;


open FIN, "gunzip -c $tgdir$tgfile | ";

my $past_region =-1;
while(($l=<FIN>)  && ($past_region==-1)) {
    $l =~ s/[\r,\n,\l]+$//g;
    if(!($l =~ /^#.*/)) {
	 my @data = split(/\t/,$l);
	 my $chr = $data[$ichr];
	 if($chr == $region_chr) {
     
	     my $pos = $data[$ipos];
	     if($pos>$last_pos) {
		 $past_region=1;
	     }
	     if(($pos>=$first_pos) && ($pos <=$last_pos)) {
		 my $snpid = $data[$isnpid];
		 my $f=$data[$iinfo];	      
		 my @info = split(";",$f);
		 my $ac=0;
		 my $an=0;
		 for(my $i=0;$i<=$#info;$i++) {
		     my @kv = split("=",$info[$i]);
		     if($kv[0] =~ /AC/) {
			 $ac=$kv[1]; # Non-reference Allele Count
		     }elsif($kv[0] =~ /AN/) { # Total Chromosome Count
			 $an=$kv[1];
		     }
		 }
		 my $maf=0;
		 if($an>0) {
			 $maf=$ac/$an;
		 }
		 my @refallele=();
		 my @altallele=();
		 my $na=0; # Number of ambigous;

		 # Process currently read SNP
		 my @genos=();
		 for(my $i = $ifirsts;$i<=$#data;$i++) {
		     my @genoarr = split(":",$data[$i]);
		     my $genostring = $genoarr[0];
		     # SAVE INFORMATION
			if($genostring =~ /^([0-2])\|([0-2])$/) {
			    my $al1 =$1;
			    my $al2 =$2;
			    push(@genos,$al1+$al2);
			    if($al1 ==0) {
				push(@refallele,1);
				push(@altallele,0);
			    } else {
				push(@refallele,0);
				push(@altallele,1);
			    }
			    if($al2 ==0) {
				push(@refallele,1);
				push(@altallele,0);
			    } else {
				push(@refallele,0);
				push(@altallele,1);
			    }
			} else { # Ambigous phase (e.g. genotype called as "/"
			    push(@refallele,0);
			    push(@altallele,0);
			    $na++;
			    if($genostring =~ /^([0-2])\/([0-2])$/) {
				push(@genos,$1+$2);
			    } else {
				push(@genos,-1); # Missing data
			    }
			}
		 }
		 # Loop over target SNPs
#		 for(my $j=0;$j<=$#region_snps;$j++) {
#		     my $target_snpid = $region_snps[$j];
		     my $target_snpid = "Haplotype";
#		     if(exists $INFO{$target_snpid}) {
#			 my $info_ref = $INFO{$target_snpid};
#			 my @tsnpinfo = @$info_ref; # ($ac,$an,$na,$data[$iref],$data[$ialt],$maf,$tpos,$tchr); 
#			 my $tpos = $tsnpinfo[6];
#			 if(abs($tpos-$pos)<=$range) {
			     # SNPs are close enough.

#			     my $ref_ref = $REFS{$target_snpid};
#			     my $alt_ref = $ALTS{$target_snpid};
#			     my $info_ref = $INFO{$target_snpid};
#			     my @trefallele = @$ref_ref;
#			     my @taltallele = @$alt_ref;
#			     my $tgenos_ref = $GENO{$target_snpid};
#			     my @tgenos = @$tgenos_ref;
			     
			     my $n_tref_ref =0; # a
			     my $n_tref_alt=0; #c
			     my $n_talt_ref=0; #b
			     my $n_talt_alt=0; #d
			     my $n_tref_na=0;
			     my $n_talt_na=0;
			     my $n_na_ref=0;
			     my $n_na_alt=0;
			     my $n_na_na=0;
			     # Loop over all the phase alleles for the samples
			     for(my $k=0;$k<=$#refallele;$k++) {
				 my $hasHapl = $hasHaplotype[$k];
				 if($refallele[$k]==1 && $hasHapl==1) {
				     $n_tref_ref++;
				 } elsif($altallele[$k]==1 && $hasHapl==1) {
				     $n_tref_alt++;
				 } elsif($refallele[$k]==1 && $hasHapl==0) {
				     $n_talt_ref++;
				 } elsif($altallele[$k]==1 && $hasHapl==0) {
				     $n_talt_alt++;
				 } elsif($hasHapl==-1) {
				     if($refallele[$k]==0 && $altallele[$k]==0) {
					 $n_na_na++;
				     } elsif ($refallele[$k]==0 && $altallele[$k]==1) {
					 $n_na_alt++;
				     } else {
					 $n_na_ref++;
				     }
				 } elsif ( $altallele[$k]==0 &&  $refallele[$k]==0) {
				     if ($hasHapl==0) {
					 $n_talt_na++;
				     } else {
					 $n_tref_na++; # cannot be $n_na_na .. since would have been caught earlier
				     }
				 }
			     } # End of Loop over Phased Alleles
			     # Two by Two Contengency table of current Allele (Vertical) vs target Alleles (new Phenotype) Horizontal
			     my $a = $n_tref_ref; # a
			     my $c = $n_tref_alt; #c
			     my $b = $n_talt_ref; #b
			     my $d = $n_talt_alt; #d
			     my $a_plus_b = $a+$b;
			     my $a_plus_c = $a+$c;
			     my $n = $a+$b+$c+$d;
			     my $chisq = ChiSquare($n_tref_ref,$n_talt_ref,$n_tref_alt,$n_talt_alt);
			     my $chisq_pvalue = ChiSquareToPvalue($chisq);
			     my $fisher_two_tailed = calculateStatistic(n11=>$n_tref_ref,n1p=>$a_plus_b,np1=>$a_plus_c,npp=>$n);
			     # Does the NA in the alt-allele dichotomize the target allele?
			     my $a = $n_tref_na;
			     my $c=  $n_talt_na;
			     my $b = $n_tref_ref + $n_tref_alt; 
			     my $d = $n_talt_ref + $n_talt_alt; 
			     my $fisher_missing_in_test = calculateStatistic(n11=>$a,n1p=>($a+$b),np1=>($a+$c),npp=>($a+$b+$c+$d));
			     
			     # Does the NA in the target allele dichotomize the test alleles.
			     
			     $a =$n_na_ref;
			     $b =$n_na_alt=0;
			     $c =$n_talt_ref+$n_tref_ref;
			     $d =$n_talt_alt+$n_tref_alt;
			     my $fisher_missing_in_target = calculateStatistic(n11=>$a,n1p=>($a+$b),np1=>($a+$c),npp=>($a+$b+$c+$d));
			     
			     # Does Missing in target lead to missing in test SNP
			     $a=$n_na_na;
			     $b=$n_na_ref+$n_na_alt;
			     $c=$n_tref_na+$n_talt_na;
			     $d=$n_tref_ref+$n_tref_ref+$n_talt_ref+$n_talt_alt;
			     my $fisher_missing_causes_missing = calculateStatistic(n11=>$a,n1p=>($a+$b),np1=>($a+$c),npp=>($a+$b+$c+$d));
			     
			     my $rsquare =rsquareLD($n_tref_ref,$n_talt_ref,$n_tref_alt,$n_talt_alt); # Can use R and Hypergeometric function to get a p-value
#			     my $pearson = rsquarePearson(\@tgenos,\@genos);
			     my $pearson = "NA";
			     my $tref = $REFAL{$target_snpid};
			     my $talt = $ALTAL{$target_snpid};
			     my $ref = $data[$iref];
			     my $alt = $data[$ialt];
			     #my @tsnpinfo = @$info_ref; # ($ac,$an,$na,$data[$iref],$data[$ialt],$maf,$tpos,$tchr); 
#		             my $tmaf=$tsnpinfo[5];
			     print "$target_snpid\t$tmaf\t$chr\t$tpos\t$tref\t$talt\t$snpid\t$maf\t$pos\t$ref\t$alt\t$rsquare\t$pearson\t$chisq\t$chisq_pvalue\t$fisher_two_tailed\t$fisher_missing_in_test\t$fisher_missing_in_target\t$fisher_missing_causes_missing\t$n_tref_ref\t$n_tref_alt\t$n_talt_ref\t$n_talt_alt\t$n_tref_na\t$n_talt_na\t$n_na_ref\t$n_na_alt\t$n_na_na\n";
			# } # End of "if" over distance between target and current SNP
		   #  } # End of if over targetSNPs for which we actually have data
		# } # end over Loop over "target" SNPs
	    } # End of "if" over outer boundaries of target region.
	 } # End of "if" over target chromosome
     } # Only process non-comment lines
   } # while statement over input VCF file.
close FIN;








sub ChiSquare() {
    my ($a,$b,$c,$d)= @_;
    # Calculate Horizontal Margins
    my $a_plus_b = $a+$b;
    my $c_plus_d = $c+$d;
    # Calculate Vertical Margins
    my $a_plus_c = $a+$c;
    my $b_plus_d = $b+$d;

    my $n = $a+$b+$c+$d;
    my $chisq  = "NA";
    if($n>0) {
	# OK, Can compute Chi-Square.
	my $e_a = ($a_plus_b/$n) * ($a_plus_c); # Simplified by ($a_plus_c/$n)*$n
	my $e_b = ($a_plus_b/$n) * ($b_plus_d); # ...
	my $e_c = ($c_plus_d/$n) * ($a_plus_c); # ...
	my $e_d = ($c_plus_d/$n) * ($b_plus_d); # ...
	# Chi-square not reliable if cell counts are too low.
	if(!($e_a<10 || $e_b <10 || $e_c<10 || $e_d<10)) {
	    $chisq = ($a-$e_a)/$e_a + ($b-$e_b)/$e_b + ($c-$e_c)/$e_c + ($d-$e_d)/$e_d ;
	}
    }
    return $chisq;
}


sub ChiSquareToPvalue() {
# One sided Chi-square Test.
#    my @chisqvals = (0.004,0.02,0.06,0.15,0.46,1.07,1.32,1.64,2.07,2.71,3.84,5.02,5.41, 6.6349,7.8794, 9.14, 10.8276,  12.1157, 13.41,  19.5114, 20.8373,22.1665,23.928,28.374,32.843,37.3249,41.8215);
#    my @pvalsvals = (0.95,0.90,0.80,0.70,0.50,0.30,0.25,0.20,0.15,0.10,0.05,0.025,0.02, 0.01,  0.005,0.0025, 0.001,    0.0005, 0.000025,0.00001,   0.000005, 0.000001,0.0000001,0.00000001,0.000000001,0.0000000001);



    my $chisq=shift @_;
    my $chisqprob = "NA";
    if($chisq eq "NA") {
	return "NA";
    } else {
	if($chisq>0) {
	   # Only precise to 5 digits.
	   $chisqprob=Statistics::Distributions::chisqrprob (1,$chisq);
       }
    }
    return $chisqprob;
}

# Compute the pearson r^2 between two vectors of genotypes encoded as 0,1,2 - skip missing data
sub rsquarePearson {
    my $geno1_ref = shift @_;
    my $geno2_ref = shift @_;
    my @geno1raw = @$geno1_ref;
    my @geno2raw = @$geno2_ref;
    my @geno1=();
    my @geno2=();
# Filter out missing data in either sample
    for(my $i=0;$i<=$#geno1raw;$i++) {
	if($geno1raw[$i]>=0 && $geno2raw[$i]>=0) {
	    push(@geno1,$geno1raw[$i]);
	    push(@geno2,$geno2raw[$i]);
	}
    }
    my $r2="NA";
    my $n=1+$#geno1;
    if($n>1) { # Need at least two samples
	my $avg1 = (sum(\@geno1))/$n;
	my $avg2 = (sum(\@geno2))/$n;
	my $cov=0;
	my $sd1=0;
	my $sd2=0;
	for(my $i=0;$i<$n;$i++) {
	    my $g1 = $geno1[$i];
	    my $g2 = $geno2[$i];
	    my $d1 = $g1-$avg1;
	    my $d2 = $g2-$avg2;
	    my $d1sq = $d1* $d1;
	    my $d2sq = $d2* $d2;
	    $sd1+=$d1sq;
	    $sd2+=$d2sq;
	    $cov+=($d1sq*$d2sq);
	}
	if($sd1>0 && $sd2>0) {
	    $r2=$cov/($sd1*$sd2);
	}
    }
    return $r2;
}

# Calculate LD r^2
sub rsquareLD {
    my $n11=shift @_;
    my $n12=shift @_;
    my $n21=shift @_;
    my $n22=shift @_;
    
    my $r2="NA";
    my $n=$n11+$n12+$n21+$n22;
    if($n>1) { # Need at least two samples
	my $det = $n11*$n22-$n12*$n21;
	my $denom = ($n11+$n12)*($n21+$n22)*($n11+$n21)*($n12+$n22);
	if($denom>0) {
	    $r2 = $det*$det/$denom;
	}
    }
    return $r2;
}


sub sum() {

    my $arr_ref = shift @_;
    my @arr = @$arr_ref;
    my $s=0;
    for(my $i=0;$i<=$#arr;$i++) {
	$s+=$arr[$i];
    }
    return $s;
}

##fileformat=VCFv3.3
##INFO=DP,1,Integer,"Total Depth"
##INFO=HM2,0,Flag,"HapMap2 membership"
##INFO=HM3,0,Flag,"HapMap3 membership"
##INFO=AA,1,String,"Ancestral Allele, ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/technical/reference/ancestral_alignments/README"
##reference=human_b36_both.fasta
##FORMAT=GT,1,String,"Genotype"
##FORMAT=DP,1,Integer,"Read Depth"
##FORMAT=CB,1,String,"Called by S(Sanger), M(UMich), B(BI)"
##rsIDs=dbSNP b129 mapped to NCBI 36.3, August 10, 2009
##INFO=AC,-1,Integer,"Allele count in genotypes"
##INFO=AN,1,Integer,"Total number of alleles in called genotypes"
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA06985	NA06986	NA06994	NA07000	NA07037	NA07051	NA07346	NA07347	NA07357	NA10847	NA10851	NA11829	NA11830	NA11831	NA11832	NA11840	NA11881	NA11894	NA11918	NA11919	NA11920	NA11931	NA11992	NA11993	NA11994	NA11995	NA12003	NA12004	NA12005	NA12006	NA12043	NA12044	NA12045	NA12144	NA12154	NA12155	NA12156	NA12234	NA12249	NA12287	NA12414	NA12489	NA12716	NA12717	NA12749	NA12750	NA12751	NA12760	NA12761	NA12762	NA12763	NA12776	NA12812	NA12813	NA12814	NA12815	NA12828	NA12872	NA12873	NA12874
#//not commented in original source file 1	533	.	G	C	-1.0	0	AA=.;AC=6;AN=120;DP=423	GT:DP:CB	0|0:6:SMB	0|0:14:SMB	0|0:4:SMB	0|0:3:SMB	0|0:7:SMB	0|0:4:SMB	0/1:6:MB	0|0:3:SMB	0|0:13:SMB	0|0:1:SMB	0|0:14:SMB	0|0:10:SMB	0|0:6:SB	0|0:2:SMB	0|0:6:SMB	0|0:4:SMB	0|0:2:SMB	0|0:15:SMB	0|0:2:SMB	0|0:1:SMB	0|0:26:SMB	0|0:6:SMB	0/1:14:MB	0|0:5:SMB	0|0:3:SMB	0|0:20:SMB	0|0:3:SMB	0|0:2:SMB	0|0:4:SMB	0|0:12:SMB	0|0:1:SMB	0|0:7:SMB	0|0:2:SMB	0|0:25:SMB	0|0:9:SMB	0/1:1:MB	0|0:9:SMB	0|0:1:SMB	0|0:6:SMB	0|0:12:SMB	0|0:7:SMB	0|0:18:SMB	0|0:2:SMB	0|0:2:SM	0|0:38:SMB	0|0:3:SM	0|0:3:SMB	0|0:5:SMB	0|0:5:SMB	0|0:3:SMB	0|0:0:MB	0|0:5:SMB	0|0:7:SMB	0|0:0:SMB	0|0:6:SMB	1|0:5:SMB	0|0:4:MB	0|0:5:SMB	0/1:5:MB	0|1:9:SMB


